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ABSTRACT 



A three dimensional primitive equation model based on 
the Boussinesq equations is developed and applied to the 
mesoscale. The model is tested with one dimensional flow in 
a balanced state and with two cases of gravity waves, which 
are forecast for 30 minutes and compared with the analytic 
solutions . 
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I. 



INTRODUCTION 



Several recent attempts to understand the complex 
dynamic interaction of the thunderstorm and its environment 
have been made using the dynamical thunderstorm model pre- 
sented by Newton and Newton (1959) . The model proposes that 
strong vertical drafts within the thunderstorm create an 
effective barrier to the environmental flow. The analysis 
by Browning (1964) concerning airflow near and within severe 
thunderstorms provides valuable additional knowledge about 
the character of the vertical motion. Convective storm mass 
and water budget assessments such as the one conducted by 
Newton (1966) substantiate the concept of upper level 
entrainment of cool, dry air first recognized by Normand 
(1946) . 

Radar determined chaff trajectories were used success- 
fully by Fankhauser (1968) to follow the motion of air 
parcels in and around a thunderstorm. Further use of chaff 
will undoubtedly make possible a good understanding of the 
paths taken by environmental parcels from all positions 
around a local storm. 

As noted by Lilly (1962) , the evolution of convective 
motions with relatively large amplitudes can best be 
described through the combined application of conventional 
analytic methods and numerical experimentation. Ogura (1963) 
discussed the possible use of the primitive equations as a 
basis for forecasting small scale dynamical features. He 
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suggested that interaction between convective elements and 
the prevailing wind field could thus be investigated numeri- 
cally. The results of such efforts, correlated with the 
subjective modeling, should produce verifying as well as 
augmenting information for the subjective products. 

A numerical study of turbulent flow conducted by 
Deardorff (1970) makes use of a three dimensional primitive 
equation model. The forecast procedure he uses is similar 
to the procedure used in this study. 

The numerical model presented here is a three dimen- 
sional primitive equation model based on the Boussinesq 
equations. Forecasts are made for the wind components and 
potential temperature after imposing initial conditions on 
these same parameters. 

The grid used for computations is 25 by 25 in the 
horizontal, with 15 levels in the vertical. Grid interval 
in all directions is 1 kilometer. A mean latitude of 35° 
is used to determine the Coriolis parameter, which is 
treated as a constant. 

Predictions are made using the IBM OS/360 computer of 
the W. R. Church computer center. A one minute time step 
was chosen to conform to computational stability require- 
ments . 
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II. DEVELOPMENT OF THE MODEL 



A. THE FORECAST EQUATIONS 

The pressure term in the equation of motion may be 
expressed as 



p v 3 p = — v 



RT 

F V 3 



1/rt 



g 



1 A< 



= • © V 3 $ , 



where 
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Using this result in the component equations of motion 
yields 
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Forming the perturbation expressions 

g = g D (z) + g' (x,y,z, t) 



and 
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© = 0 O + 0 (x,y, z, t) 



where 



I $'\ « Bo and [ 0 | « 0 o ' 



and substituting into (1) yields 



du 
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where the relatively small product 
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is neglected. 
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For the vertical component equation the substitution yields 
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0 are related by the hydrostatic equation 
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the vertical component equation may be written 
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Defining <p = Q q ft' and substituting into (4), (5) and 

(6) results in the system of equations 
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continuity equation. In order to eliminate sound waves, 
incompressible flow is assumed. Ogura and Phillips (1962) 
have shown that this system of equations, which they call 
the anelastic equations, is quite accurate for shallow con- 
vection. 

B. DETERMINATION OF THE PRESSURE TERMS 

Equations (7) , (8) , and (9) may be combined and written 

in the vector form 



a t 



= -(V V V 3 )V, - Vo<*> -f (k X Vo) + -T— k + v V 
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(12) 



Taking the three dimensional divergence of (12) yields 
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For a limited horizontal scale the Coriolis parameter may be 
treated as constant. Thus, (13) becomes 



V 3 2 * 




(k x v 3 ) 
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(14) 



which can be solved using relaxation techniques. 
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III. FINITE DIFFERENCING 



Equations (7) through (10) are solved numerically by 
introducing finite differences in x, y, z and t. Solution 
for the pressure term can be carried out by applying the 
extrapolated Liebmann relaxation technique to evaluate (14) . 
Appendix A shows the development of the relaxation equation. 
The procedure used to compute the forcing function is 
included in Appendix B. 

In order to maintain finite differencing consistency 
throughout the entire forecast process, it is necessary to 
use the finite difference forms of equations (7), (8) and 

Q 

(9) to form the finite difference equivalent of equation 
(14) . 

The vertical domain consists of fifteen levels with 
values of u, v, w, Q, and <p assigned to each grid point. 

The space differencing method used is based on a scheme 
devised by Arakawa (1966) and is designed to conserve total 
energy. The non-linear operator, £(s) , takes the finite 
difference form 



(S) “ [ (u i + l,j.k + + 



+ 



[^ V i,j+l/h + v i,j,k^ ^ S i/ 

(vi,j-l,k + v i/j/k ) (Si,j 



j+l,k + s i, j,k } 

-l,k + s i ( j |k )] /4A ^ 
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[ (w i, j,k+l 



+ w i,j,k^ ^i,j,k+l + 



(w i,j,k-l + w i, j , k } (S i,j,k-l + S i < j,k ) ] /4 A: 



The linear space derivatives are of the centered form 

as = (S i+l,j,k ~ S i-l,j,k) 

dx 2 Ax 

Centered time differencing is used for all forecasts except 
the first which employs a forward time step. 
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IV. BOUNDARY CONDITIONS 



The boundary conditions are selected to enforce mass and 
energy conservation within the domain of the grid. Figure 1 
illustrates the grid domain and boundary placement in the 
horizontal. Boundaries in the vertical are similarly 
placed, with bottom and top boundaries at 



k 



2+3 

2 



2.5 



and 



k = 



(KM-1) + (KM-2) 
2 



respectively, where KM is the number of grid points in the 
vertical. 



i=l 

j=JM 



i=IM 

j=JM 



» ' « 

I 



North boundary 






West! boundary 



East [boundary 



• • 



1 



South boundary 



i=l 

j=l 

Figure 1: Horizontal grid 



i=IM 

j=l 



Periodicity in the east-west direction ensures cyclic 
continuity of all parameters. This condition is imposed by 
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applying the relationships 




S IM-3,j,k ' 



S 



2, j ,k 



S 



IM-2,j,k ' 




and 



S IM,j,k ^4 , j , k 



where & = u, v, w, 0 and <p . Since values of S at points 
in the outer columns are required in solving for the pressure 
term and in calculating divergence, a four point periodicity 
constraint is necessary. 

In order to prevent the flow of mass and energy through 
the northern and southern boundaries, a condition of no flux 
is imposed on the v component of the wind. This is accom- 
plished by 



v i, 2 , k " " v i,3,k 



and 



v . 



i, JM-l,k 



v i , JM-2 , k 



The values of v on the outer row are obtained by 
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The u and w wind components and 0 are projected 
outward by 



S , — S . on — S . on 

i/l/X i / 2 # k i/ 3 / k 



and 



S i, JM-1, k S i,JM-2,k ' 



where S = u, w, 0 



Values of pressure on the outer rows are obtained by using 
the geostrophic relationship 



%<t> 



= -fu 



The upper and lower boundaries are treated in the same 
manner as the north and south boundaries, with the no-flux 
condition being imposed on w by 



w. , 0 = -Wj; i o , 

l, j , 2 J ' J 



W i,j,KM-l _w i, j ,KM-2 , 



W. . n = W. . . , 

l / J / 1 1/ J / 4 



and 



W i,j,KM W i , j , KM- 3 



Values of u and v on the top and bottom levels are 
obtained by 

S i,j,l = S i , j , 2 = S i , j , 3 
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and 



S i,j,KM “ S i/ j , KM-1 “ S i,j,KM-2 

where s = u, v. 

The hydrostatic relationship, 

b<P = 2®. 

Sz Oo 

is used to determine values of 0 and 0 above the top and 
below the bottom boundaries. The application of proper 
boundary conditions was found to be critical. Appendices C 
and D show the development of these conditions in detail. 
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V. INITIAL CONDITIONS 



Initialization of the u,v,w and 0 fields was 
accomplished according to the experiment being prepared. 
During the check-out phase of programming the model, a 
steady-state, geostrophic, zonal flow was used to provide a 
check on computational stability. Subsequent tests were 
initialized to include gravity waves, both stable and 
unstable . 

A. ZONAL FLOW CASE 

Stability of the balanced state was tested by inserting 
a 10 m sec~^ wind throughout the u field. Initial v 
and w fields were set to zero. Standard atmospheric 
values of pressure were introduced and, by using an adia- 
batic lapse rate of temperature, the equation 

r 1 k 




provided the initial 0 field. Since the relaxation 
technique used to obtain the initial <t> field proved to be 
extremely efficient, choosing the initial guess <p field as 
an arbitrary constant was convenient. 

B. GRAVITY WAVE CASES 

Linearized gravity wave solutions, accurate for small 

y 

scale perturbations only, were used to test the model. For 
the stable gravity wave the following relationships describe 
the initial conditions. 
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u = — 



A m w 

A 2 + /r 



sin Ax cos 



UY 



cos mz 
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fl m W 



cos Ax sin VY cos 



mz 



w 



= W cos Ax cos IXy sin mz 



W d 0 

Ac £ z 



in Ax cos ^(y sin mz 



sm 



Phase speed may be calculated from the relationship 



= h J A 2 + fJ 
A » A 2 + H 2 + m 2 



( 15 ) 



For the gravity wave associated with unstable stratifi- 
cation, the initializing relationships are 
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and 



-N 2 (A 2 + I U 2 ) 

X 2 + H 2 + m 2 



( 16 ) 



where n is the growth rate of the disturbance. 

The <f> field may be initialized for the gravity waves 
exactly as in the zonal flow case. 
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VI. RESULTS 



The balanced state introduced for the one dimensional 
flow case as a check on the computational stability and 
forecast capability of the model should remain indefinitely, 
subject to the validity of the approximations in the basic 
equations, round-off error in the computer, accuracy of the 
assumptions at the boundaries, and stability requirements. 
The use of proper boundary conditions was the most important 
factor in determining the length of the forecast produced. 
Until the combination of boundary conditions shown in 
section IV and Appendix C was utilized, the forecasts were 
valid for no longer than 5 minutes. Longest forecasts of 
the balanced flow were obtained by applying the boundary 
conditions that are periodic in x, geostrophic in y, and 
hydrostatic in z. 

Accuracy of forecasts was monitored through the calcula- 
tion of the three dimensional divergence, which should have 
remained very small throughout the prognosis. The smallest 
possible relaxation tolerance consistent with single preci- 
sion mode capability of the computer was found to allow the 
least amount of divergence. When larger tolerances were 
allowed, a substantial amount of divergence developed on the 
north and south boundaries at mid-levels and increased with 
he ight. 

A test was conducted to determine the optimum relaxation 
coefficient for the grid interval and grid domain which the 
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model employs. The value selected did not appear too 
critical, as long as it was between 1.5 and 1.7. 

For the first experiment that produced a relatively long 
prognosis, an adiabatic lapse rate of temperature and 
standard atmospheric values of pressure were used to ini- 
tialize the potential temperature field. This produced a 
slightly unstable average lapse rate of potential tempera- 
ture, with two thin stable layers included near the bottom. 
With these conditions the model forecast for 47 minutes. 

The time step was then decreased from 1 minute to 40 
seconds, with the result that the forecast proceeded for 72 
time steps, which was essentially the same length as the 
previous forecast. In an effort to improve the initial 
balance, pressure values which would lead to a nearly neutral 
lapse rate of potential temperature were inserted in the 
first experiment. The forecast continued for 56 minutes. 
Returning again to the conditions used in experiment one, 
the lapse rate of potential temperature was modified to 
include stable layers across the top and bottom boundaries. 
With these conditions the model produced a forecast for 69 
minutes . 

Introduction of a small perturbation in the form of a 
gravity v/ave provided another method of observing the model 
behavior. The gravity wave in stable stratification should 
propagate in the x direction only, without growing. With 
unstable stratification the wave should not propagate, and 
should grow at a rate given by equation (16) . 
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The wave inserted for the stable case was forecast for 
30 minutes, after which a comparison was made between the 
initial and forecast values of u, v, w and 0 . With this 
comparison it was possible to observe the amount of spurious 
growth and the speed of propagation. The relationship of 
the initial and forecast wave forms in x for the line j=8, 
k=6 is shown in Figures 2a, 2b, 2c and 2d. The exact solu- 
tion of equation (15) indicated that the phase speed for the 
case shown was 10.2 m sec - "*". Phase speeds calculated from 
the plotted values were 10.5 m sec"'*', 10.4 m sec”-*- , 

9.5 m sec ^ and 10.3 m sec -1 for u, v r w and 0 respec- 
tively. The difference in phase speeds between w and 9 
was in no doubt related to the fact that the wave grew. The 
wave forms in y illustrated in Figures 3a, 3b, 3c and 3d 
show that the wave grew in the y direction but did not 
propagate. Figures 4a, 4b, 4c and 4d show that no propaga- 
tion of the wave occurred in the z direction, but that the 
growth was not smooth. 

For each combination of initial conditions applied to 
the wave, growth was apparent. By modifying the lapse rate 
of potential temperature toward adiabatic conditions, the 
growth decreased. 

The gravity wave with unstable stratification was 
compared in the same manner as was the stable case. Figures 
5a, 5b, 5c and 5d show that the wave did not propagate. 
Furthermore, the growth that was experienced by the wave was 
very close to that which was expected. The amplification 
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factor based on the computed growth rate of .1376 x 10 - ^ 
sec was determined to be 11.9 at 30 minutes, and amplifi- 
cation factors of 11.8, 12o0, 12.4 and 11.3 were noted for 
u, v, v/ and 6 respectively. 
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Fig. 2a — The u component of the wave form in the x direction for 
k=6, stable case 
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Fig. 2b - The v component of the wave form in the x direction for 
k=6, stable case 
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Fig. 2d - The 0 wave form in the x direction for j=8, k=6, stable case 
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Fig. 3a - The u component of the wave form in the y direction 
stable case 
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Fig. 3c - The w component of the wave form in the y direction for 
k=6, c table case 
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Fig. 3d - The 0 wave form in the -- direction for i=6, k=G, stable ca 
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Fig. 4a - The u component of the wave form in the z direction for 
i=8, stable case 
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Fig. 4b - The v component of the wave form in the z direction for 
i=8, stable case 



o 

ro 

II 

-P 





X 

0) 

"d 

a 

•H 



— 






in 



— ro 



o 

00 



o o 



T 

o 

i 



T 

o 

CO 

I 



( fr _0T x T- oas ai) m 



40 



Fig. 4c - The w component of the wave form in the z direction for 
i=6 , i=8, stable case 
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Fig. 4d - The 9 wave form in the z direction for 
i=6, j=8, stable case 
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Figo 5a - The u component of the wave form in the x 
direction for j=8, k=6, unstable case 
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Fig. 5b - The v component of the wave form in the x 
direction for j=8, k=6, unstable case 
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- The w component of the wave form in the x 
direction for j=8, k=6, unstable case 
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VII. SUMMARY AND CONCLUSIONS 



A multi-layer primitive equation model was applied to 
a mesoscale grid for two test situations. First, a one 
dimensional, balanced, steady state current was initialized 
to check the computational stability and finite difference 
form of the equations of motion 0 Then an initial gravity 
wave perturbation was introduced for both stable and 
unstable stratification. 

The stable stratification case yielded a reasonably 
accurate phase speed of propagation but some undesirable 
growth of the disturbance was observed, as was discussed in 
the previous section. 

The unstable stratification case produced nearly correct 
growth rates but some difference in the w and d growth 
rates may be an indication that a problem exists similar to 
that in the stable case, where w and 6 waves were propa- 
gated at slightly different speeds. Apparently, some small 
error still exists in the program. 

The fact that the growth of the propagating gravity 
wave can be decreased by inserting a more nearly adiabatic 
lapse rate of potential temperature can be explained by 
observing that the buoyancy term in the Boussinesq equations 
becomes increasingly inaccurate with increased height. By 
assigning to 6> 0 a value of potential temperature from a 
middle level, the error introduced in the buoyancy term is 
minimized. 
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Since the Boussinesq equations are accurate only for a 
shallow atmosphere, the model could be expected to forecast 
more accurately if the depth were decreased. In such a case 
the problem with the buoyancy term would be further mini- 
mized. 

Further studies using this model should have a provi- 
sion for including vertical compressibility. Additionally, 
since a primary source of the energy in a thunder storm- type 
cloud is latent heat, application to cumulus scale convec- 
tive activity would be incomplete without the inclusion of 
moisture„ 
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APPENDIX A 



EXTRAPOLATED LIEBMANH RELAXATION SCHEME FOR 
THREE DIMENSIONAL APPLICATION 



The equation to be solved by relaxation is 



V 3 <P - F = 0 



(A-l) 



where 



F = 



[• 



- V I (V 3 ‘ V 3 )V 3 I + £ ° f + 



] 



g bd 

6q d z 



Writing (A-l) in finite difference form and expanding yields 



<t>. ± 0 . , - 20. . . + Cf>. n , 

i+2,j,k ^ i, j,k ^ ^1-2, j,k 

(2 A x) 2 



ftj, j + 2,k ~ 2 4>1, j,k + ^i, j -2 , k 

(2 Ay) 2 

+ 1 /k+2 ~ 2 J ' k + ^i/ j , k-2 

(2 Az) 2 



- F. . , = 0 
1/ J /h 



Rearranging terms leads to 

(2Ay)2 (2 Az)2 <‘f > i+2,j,k + 

+ (2 Ax) 2 (2 Az) 2 (<*>i /j + 2/ k + <*>1^-2,^ 



+ (2 Ax) 2 (2 Ay) 2 (<*> i/j/k+2 + < ^ > i/ j / k-2 ) 
- (2 Ax) 2 (2 Ay) 2 (2 Az) 2 F if j <k J A 
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where A is the fraction 



2 (2Ay) 2 (2Az) 2 + (2 Ax) 2 (2Az) 2 + (2Ax) 2 (2Ay) 2 



1 



old 



Adding and subtracting <t> produces 

i, j ,k 



new 



old 



„ = <>> . . , + A [ (2 Ay) 2 (2 A z) 2 ( <t>, . + <f> ) 

i ,j,k i,j,k L i - z , j,k 



+ (2Ax) 2 <2Az> 2 < <*>._. +2>k + j. 2jk ) 



+ (2 Ax) 2 (2Ay) 2 < <?>._. (k+2 ^ i(j , k _ 2 ) 

$ 0ld 

- (2 Ax) 2 (2 Ay) 2 (2Az) 2 F i( j ( J 

A J 



For successive over relaxation the coefficient Cl is 
included, which results in the final form 



new 

0 . - (i 

i,j.,h U 



old 

-0) + a A [ (2Ay) 2 (2Az) 2 ( <0 i+2ij/ 



+ ^i-2,j,k^ 



+ < 2A x) 2 (2A 2 ) 2 (0 i/j+2 , k + 0 i, j - 2 ,k) 

+ (2 Ax) 2 (2Ay) 2 ( 0 i/ j / k+2 + 0i,j ,k-2> 

- (2 Ax) 2 (2 Ay) 2 (2 Az) 2 (A-2) 
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Note that if Ax = Ay = Az, equation (A-2) reduces to the 
more familiar 



new 



old 



^i/j/h ^ ^i/j/k + 6 [ ^i+2 , j , k + ^i-2,j,k 



+ < ^ > i / j+2,k + ^ > i,j-2,k + ^^^+2 + < ^ > i,j f k-2 



- (2 Ax) F 



i#j/k] * 
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APPENDIX B 



DEVELOPMENT OF THE FORCING FUNCTION 



In order to maintain consistent finite differencing 
between the prognostic equations and the forcing function of 
the Liebmann relaxation equation, a procedure which was 
recommended by Dr. R 0 T. Williams (personal communication) 
is used to develop the forcing function. The process is 
carried out in 6 steps for each point in the grid. The 
6 steps required to build F 6 ^ 5 the forcing function 

at i=5, j=5, k=5 , are 



,1 , ^^ 5 , 5,4 

5,5,5 2 Az 



2 — G ( v) 545 

+ F 1 

5,5,5 2 Ay 5,5,5 



3 _ - g < u >4,5,5 2 

5,5,5 2 A x 5,5,5 



p 4 , G(u) 6,S,5 3 

5,5,5 2Ax 5,5,5 



5 g < v >5,6,5 4 

F = 1 — L — + F 

5,5,5 2 Ay 5,5,5 



r G (w) 555 

? b - + f 5 

5,5,5 2 Az 5,5,5 



where 
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and 



G (w) = - 


£ (w) 


G (v) = - 


£(v) 


G (u) = - 


£ (u) 



4 - ° 2 

+ % + * V w, 

" fu + V V 3 2 v, 

+ f v + y u. 
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APPENDIX C 



BOUNDARY CONDITIONS AT THE WALLS 

For the plane described by j=3, the equations of motion 
may be written 



a u 

St 

i/ 3,k 



i, 3,k 



S0 

* x i,3,k 




3 , k 



(C-l) 



S v 



S t 

i z 3 / k 



£ (v) 



i/ 3 , k 



S0 

Sy 

i, 3 , k 



■^ u i, 3,k 



(C-2) 



Sw 
S t 



i/ 3,k 



£ (w) 



i, 3,k 



Sz 



+ 

i, 3,k 




3,k 




(C — 3 ) 



Forming the time derivative of divergence from (C-l ) , (C-2) , 

and (C-3) and using finite differencing yields 



0 



S u _ S u 

d t i+ 1 , 3 , k i-i, 3 1 k 

2 Ax 



S v _ S v 

i, 4 , k i/ 2,k 

T 

2 A y 

S w £ w 

at i,j,k+l St i / j / k-1 

+ 

2 A z 



(C-4) 
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mtroaucin? 



the 



con 






v i, 3 ,* 






„*ic* » lieS 



dj£ 

n i( 2,^ 



J^v. 

^ i,3,* 





Applying the conditions 



u i,2,k “ u i,3,k 



(C-8) 




(C-9) 



and 




(C-10) 



produces 





which is the equation to be solved by relaxation. 

Thus, the conditions for the southern boundary must 
include the relationships given by (C-5) , (C-8) , (C-9) and 

(C-10) . 

A similar approach produces corresponding conditions 
at the northern boundary. These conditions are 



v i, JM-l,k “ “ v i,JM-2,k 



(C-ll) 



u i,JM-l,k ~ u i,jM-2,k 



(C— 12) 



0 . , - 0 . 

i, JM, k ^i,JM-2,k 



-2 Ayf u- 



i, JM-l,k 



(C— 13) 



and 




-2 Ayf u 
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For the plane described by k=3, the same process can 
be carried out to develop boundary conditions at the top and 
bottom, which are 



w - -• - “ -w i , j , 3 ' 



±, 3,2 



(C-15) 



<?> . - d>. = a 

i/j/3 i,j,l 0 o ±, 3,2 ' 



(C-16) 



<t> . <b . = 2 .4 Z g. „ 

i/j/4 i,j,2 0 ^ i» j / 



(C-17 ) 



w- 



i, j , KM— 1 “ -w i , j , KM- 2 ' 



(C-18) 



<t>. . 

i, J/KM 






j , KM- 2 



2 Azg 



^ i, j , KM-1 



(C-19) 



and 



<t> 

i, j , KM-1 



<t> 



i , j , KM-3 



2 . 4 . z g, 9 

0 O / j / KM— 2 * 

(C-20) 
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